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ABSTRACT 

(— I \ We present an analysis of the effects of baryon physics on the halo mass function. 

The analysis is based on simulations of a cosmological volume having a comoving size 
of 410 /i-i Mpc, which have been carried out with the Tree-PM/SPH GADGET-3 
code, for a WMAP-7 ACDM cosmological model. Besides a Dark Matter (DM) only 
simulation, we also carry out two hydrodynamical simulations: the first one includes 
^ I non-radiative physics, with gas heated only by gravitational processes; the second 

one includes radiative cooling, star formation and kinetic feedback in the form of 
galactic ejecta triggered by supernova explosions. All simulations follow the evolution 
of two populations of 1024'^ particles each, with mass ratio such to reproduce the 
assumed baryon density parameter, with the population of lighter particles assumed 
to be coUisional in the hydrodynamical runs. We identified halos using a spherical 
overdensity algorithm and their masses are computed at three different overdensities 

■ (with respect to the critical one), Ac — 200, 500 and 1500. 

We find the fractional difference between halo masses in the hydrodynamical and 

in the DM simulations to be almost constant, at least for halos more massive than 
log(MA^//i~^ Mq) ^ 13.5. In this range, mass increase in the hydrodynamical sim- 
ulations is of about 4-5 per cent at Ac = 500 and ~ 1 - 2 per cent at Ac = 200. 
Quite interestingly, these differences are nearly the same for both radiative and non- 
radiative simulations. Mass variations depends on halo mass and physics included for 
, higher overdensity, Ac — 1500, and smaller masses. Such variations of halo masses in- 

■ duce corresponding variations of the halo mass function (HMF). At z = 0, the HMFs 

■ for GH and CSF simulations are close to the DM one, with differences of ^ 3 per cent 

at Ac = 200, and ~ 7 per cent at Ac = 500, with ~ 10 - 20 per cent differences 
reached at Ac — 1500. At this higher overdensity, the increase of the HMF for the 
radiative case is larger by about a factor 2 with respect to the non-radiative case. As- 
suming a constant mass shift to rescale the HMF from the hydrodynamic to the DM 
simulations, brings the HMF difference with respect to the DM case to be consistent 
with zero, with a scatter of ^ 3 per cent at Ac = 500 and ^ 2 per cent at Ac = 200. 

Our results have interesting implications to assess uncertainties in the mass func- 
tion calibration associated to the uncertain baryon physics, in view of cosmological 
applications of future large surveys of galaxy clusters. 

Key words: clusters: cosmology: theory - dark matter - galaxies: formation - halos 
- methods: numerical 



1 INTRODUCTION 



An accurate calibration of the halo mass function is at the 
* wgcui@oats.inaf.it hearth of a range of cosmic structure formation studies, from 
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the stud y of g alaxy formation through semi-analytical mod- 
els (e.g. iBaugh 2006), to the co smological application of 
galaxy clusters (| Allen et al.ll201ll '). Under the standard hi- 
erarchical ACDM model, halos are formed from initial den- 
sity peaks through gravitational instability. The halo mass 
function (HMF hereafter) is directly connected to the pri- 
mordial density field. Since the abundance of density peaks 
over a given mass scale A/ only depends on the r.m.s. value 
aM of the linear fluctuation field at that mass scale, the 
abundance of halos is expected to be universal once ex- 
pressed as a function of ctm, as assumed by the Press- 
Schechter approach base d on the spherical collapse model 
JPress fc Schechtedl 19741) an d in th e ellipsoidal collapse ex- 
tension bv lSheth fc TormenI j 19991 ). 

Through the years, N-body simulations of large cosmo- 
logical volumes have been used to calibrate fltting functions 
for a universal HMF (e.g. I Jenkins et al. ll200ll : IWarre"n" et al.l 
l2006l : ISpringer et al. l l2005l : lLukic et al.ll2007l ). Thanks to the 
progressive increase in the covered dynamic range of halo 
masses, simulation results have been shown to predict sub- 
tle but sizable dev iations from unive rsality of the mass func- 
tion. For instance iReed et al.l (|2003|) found th at the univer- 
sal mass function by ISheth fc TormenI l| 19991 ) over-predicts 
the number of most massive halos found at z > 10. This re- 
sult w as confirmed by the subsequent analysis bv lReed et al.l 
iiooj), who pointed out that an even better fit for the mass 
function can be obtained if it is allowed to depend not only 
on the linear r.m.s. overdensity, but also on the local slope 
of the linear power spectrum at the relevant mass scale. Us- 
ing the spheri cal over-densi t y (SO ) algorithm to measure 
cluster masses, [Tinker et al.l (|2008h combined different sim- 
ulations to calibrate the HMF, with masses measured at dif- 
ferent overdensities. They found significant deviations from 
non-universality, with a monotonic decrease of halo abun- 
dance with increasing redshift, and provided fltting func- 
tions to such deviation. Besides confirm ing the non-universa l 
behaviour of the high end of the HMF. ICrocce et all (|2010h : 
[Tinker et al.l (2008) also pointed out that using more ac- 
curate second-order Lagrangian perturbation theory to set 
initial co nditions could is relevant fo r an accurate HMF cal- 
ibration. iBhattacharva et al.l l|201 j ) analyzed the HMF for 
an extended suite of simulations also including quintessence 
models with w —1 for the Dark Energy equation of state, 
and confirmed violation of universality at the ~ 10 per cent 
level for the range of masses and redshift covered by their 
simulations. 

At least in principle, calibrating the mass function of 
DM halos with great accuracy is just a technical problem to 
be tackled by extending the dynamic range of simulations 
and the parameter space of considered cosmological models. 
However, the back-reaction effects of baryons on dark mat- 
ter halos are known to impact on density profiles and, there- 
fore, on their mass. In turn, these back-reaction effects are 
expected to depend on the detail of the physical processes, 
such as radiative cooling, star formation and energy feedback 
from astrophysical sources, w hich determine the di stribution 
of baryons within DM halos. [Tinker et all ||2008D included 
a non-radiative hydrodynamical simulation of a large cos- 
mological volume within the large set of simulations that 
they analyzed, without however discussing in de tail the ef- 
fect of baryons on the HMF. iRudd et al.i (|2008l ) compared 
the HMF computed for a DM-only simulation with those ob- 



tained from the corresponding hydro-dynamical simulations, 
carried out with an Adaptive Mesh Refinement (AMR) code 
both with non-radiative physics and including the effect of 
gas cooling and star formation. After computing masses at 
the viral radius, they found that the HMF for non-radiative 
simulation is very close to the DM-only one, at least in the 
mass range numerically resolved by both simulations. On the 
other hand, the radiative simulation was found to produce 
a ~ 10 per cent higher mass function, as a consequence of 
the higher concentration h alo concentration res ulting from 
adiabatic contraction (e.g. iGnedin et al.l [2004) . A signifi- 
cant increase of halo concentration from adiabatic contrac- 
tion is a well known consequence of (over) efficient gas cool- 
ing (c.g.'P edrosa et al.ll2009l : lTissera et al.ll2010l : lDuffv et al.1 
2010). In line with this result on halo concentration, also 
the total matter power spectrum in radiative hydrodynamic 
simulations has been shown to have a higher amplitude 
than for DM-only N- body simulations, small non-linear 
scales k > IhM pc'^ (|Rudd et al.' '200^: I jjng et al.1 l2006l : 

~~ini;l( 



van Daalen et ahl l201ll ; ICasari ni et al.. ,20l3, ; Casarini et 
al., in preparation). However, also the simple case of non- 
radiative hydrodynamics has been suggested to increase halo 
concentration, as a consequence of a redistributi on of energy 
betw e en baryons and DM during halo collapse (iRas ia et al.l 
l2004l : iLin et al.l |2006| ). An increase of halo concentration 
turns into an incre ase of halo ma s ses, h ence increasing the 
halo mass function. IZentner et al.l (|2008l ) suggested that the 
main effects of baryons can be translated into a simple 
change of halo concentration s, thereby resu lting in a uniform 
relative shift of halo masses. IStanek et al.l ((2009, ) compared 
the cluster masses and mass functions for a set of simula- 
tions including only DM, non radiative hydrodynamics, as 
well as radiative runs with and without pre-heating. They 
reported for the pre-heated run an average decrease of halo 
mass AfsocEl by 15 percent with respect to the non-radiative 
case, and 16 percent halo mass enhancement for simulation 
with cooling and star formation (CSF) with respect to the 
DM simulation. These mass vari ations turn into diffe rences 
of the HMF of up to ~ 30 percent. [Stanek et al.l (|2009l ) based 
their analysis on two different sets of simulations, based on 
SPH and AMR codes, also using slightly different choices for 
the cosmological parameters. Furthermore, results for their 
CSF case were based only on re-simulations of the 13 most 
massive halos identified in the original simulation volume. 

In order to improve with respect to the current under- 
standing of baryon effects on the HMF, we present in this 
paper the analysis of three cosmological simulations based 
on DM only, non-radiative hydrodynamics and cooling, star 
formation and supernova (SN) feedback. These simulations 
are carried out starting from the same initial condi tions and 
using t he same Tree-PM/SPH code GADGET-3 ISpringell 
(2005a). Resolution and box-size of our simulations are ad- 
equate to cover the halo mass distribution over the range 
log(M2oo//i"^ M0) ~ (12.5 ~ 15) at z = 0. Due to the inclu- 
sion of hydrodynamics, the dynamic range covered by our 
simulations is in general narrower than that accessible by 



^ In the following, we will use the convention Ha^ to indicate the 
halo radius encompassing an average overdensity of Ac times the 
critical cosmic density pcr(^). Accordingly, is the halo mass 

contained within R/\a- 
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N-body simulations used over the last few years for preci- 
sion calibrations of the HMF. For this reason, the aim of this 
paper is not that of providing one more of such calibrations, 
rather our goal is to asses in detail the impact of baryons on 
the HMF. 

This paper is organized as follows. In Section O we de- 
scribe the simulations. Section [3] is devoted to the presen- 
tation of the analysis method and results. After describing 
the halo identification method based on spherical overden- 
sity, we present the results of our analysis in terms of mass 
variation of halos and resulting effect on the HMF. Finally, 
we discuss our results and present the main conclusions in 
Section U 



2 THE SIMULATIONS 

We carry out simulations of a flat ACDM cosmology with 
Qm = 0.24 for the matter density parameter, Qb = 0.0413 
for the baryon contribution, as = 0.8 for the power spec- 
trum normalization, — 0.96 for the primordial spectral 
index, and h — 0.73 for the Hubble parameter in units of 
100 kms~^Mpc~^. Initial conditions have been generated at 
z = 49 using the Zeldovich Approximation for a periodic 
cosmological box with comoving size L = 410 Mpc. 
Initial density and velocity fields are sampled by displac- 
ing, at redshift 2: = 41, the positions of two sets of 1024"^ 
particles each, according to the Zeldovich approximation, 
from unperturbed positions located onto two regular grids 
which are shifted by half grid size with respect to each other. 
Masses of the particles belonging to the two sets have ra- 
tio such that to reproduce the cosmic baryon fraction, with 
mi ~ 3.54 X lO^h-^ Mq and m2 ~ 7.36 x lO^h'^ Mq. In 
the DM-only simulation both particle species are treated as 
coUisionless, while in the hydro-dynamical simulations m2 
provides the mass of gas particles. We emphasize that this 
prescription to set initial conditions for the DM simulation 
ensures that it starts exactly from the same sampling of den- 
sity and velocity field as its hydro-dynamical counterpart. 
Convergence of the mass function against changing initial 
redshift and effect of using second-order Lagra ngian Pertur- 
bati on Theor y (2LPT) have bee n discussed bv I Tinker et al.l 
liobs ) and Crocce et al.1 (|2010l '). Although small but size- 
able effects have been detected in the high-end of the mass 
function, the general result is that the effect of 2LPT is 
rather small for initial redshift and resolution relevant for 
our simulations. Furthermore, since our analysis is focussed 
on the relative effect induced by the presence of baryons, we 
expect the main conclusions not to be affected by increas- 
ing the accuracy in the computation of displacements in the 
generation of initial conditions. 

Simulations are carried out using the TreePM-SPH code 
GADGET-3, an improved version of the GADGET-2 code 
l|Springel|[2005bl ). In GADGET-3 domain decomposition is 
performed by allowing disjointed segments of the Peano- 
Hilbert curve to be assigned to the same computing unit, 
thus turning into a significant improvement of the work- 
load balance when run over a large number of processors. 
Gravitational forces have been computed using a Plummer- 
equivalent softening which is fixed to epi — 7.5h~^ physical 
kpc from z — to z = 2, and fixed in comoving units at 
higher redshift. 



Besides a DM-only simulation (DM hereafter), we 
also carried out two hydrodynamical simulations. A non- 
radiative simulation only including gravitational heating of 
the gas (GH hereafter) used 64 neighbours for the computa- 
tion of hydrodynamic forces, with the width of the B-spline 
smoothing kernel allowed to reach a minimum value equal to 
half of the gravitational softening. A second hydrodynami- 
cal simulation has been carried out by including the effect of 
cooling and star formation (CSF hereafter). In this simula- 
tion radiative cool ing is computed for non-van ishing metal- 
licity according to Sutherland fc Dopital l|l993l l, also includ- 
ing heating/cooling from a spatially uniform and evolving 
UV background . Gas particles above a given threshold 
density are treated as multi-phase, so as to provide a sub- 
resolution description of t he interstellar medium, according 
to the model described bv lSpringel fc Hernguistl l|2003l ). In 
each multi-phase gas particle, a cold and a hot-phase coexist 
in pressure equilibrium, with the cold phase providing the 
reservoir of star formation. Conversion of coUisional gas par- 
ticles into coUisionless star particles proceeds in a stochastic 
way, with gas particles spawning a maximum of two gener- 
ations of star particles. The CSF simulation also includes a 
description of metal production from chemical enrichment 
contribute d by SN-II, SN -Ia and AGB stars, as described 
by (Torn<jtorrelai][2003). Kinetic feedback is implemented 
by mimicking galactic ejecta powered by SN explosions. In 
these runs, galactic winds have a mass upload proportional 
to the local star-formation rate. We use = 500 km 
for the wind velocity, which corresponds to assuming about 
unity efficiency for the conversion of energy released by SN- 
II into kinetic energy for a Salpeter IMF. The feedback 
model included in the CSF simulation is known not to be 
able to reg ulate overcooling, esp ecially in large cluster-sized 
halos (e.g. iBorgani et al.l [20041 ) . To show this, we imple- 
ment a consistent comparison in Figure [U between obser- 
vation al results on the mas s fraction in stars within R50Q 
(from ^Gonzalez et ajj l2007l : see also [Gonzalez et al.l 12007 : 
Lagana et al. 2011) and results obtained from the analysis 
of the clusters and groups identified in the CSF simulation. 
Quite apparently, simulations predict a decline of the stellar 
mass fraction as a function of cluster mass which is much 
milder that the observed one. As a result, massive systems in 
simulations are predicted to have an exceedingly high mass 
fraction in stars. Therefore, while none of the two hydrody- 
namical simulations provides a fully correct description of 
the evolution of baryons within DM halos, considering both 
the GH and the CSF runs one should provide a useful indica- 
tion of the impact of current uncertainties in the description 
of baryon physics. 



3 RESULTS 

3.1 Halo identification 

The two most common methods for halo identifications sim- 
ulations are th e one bases on th e Friend-of- Friend (FoF) 
algorithm (e.g. I Davis et al.l Il985l ) an d that based on th e 
spherical overdensity (SO) algorithm l|Lacev fc Cole|[l994l ). 
The FoF halo finder has only one parameter, b, which de- 
fines the linking length as bl where I = is the mean 
inter-particle separation, with n the mean particle number 
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Figure 1. Comparison between observational results and simu- 
lations on /», defined as the stellar mass fraction within -R50O1 
for groups and clusters of galaxies. Circles with errorbars are 
the observational resu lts for the clusters and groups analysed by 
I Gonzalez et al. 1 ||2007| '). with the dashed line showing the best-fit 
linear regression to these data points. The continous line show 
the results for the halos identified in the CSF simulations. Error- 
bars in this case refer to the r.m.s. scatter within each interval in 
M500. 



density. In the SO algorithm there is also only one free pa- 
rameter, namely the mean density Ac pcrit contained within 
the sphere within which halo mass is computed, with pcrit 
being the critical cosmic density. Each of the two halo finders 
h as its own advantages and shortcomings (see more d etails 
in I Jenkins erallbOOll : Iwhitell200ll : [Tinker et al.ll2008l . etc), 
and the difference of halo mass and HMF defined by the 
two methods have been discussed in several analysis , (e.g . 
Whit e 2002; Rccd ct al. 2003, 2007; Cohn & Whit3 1200^ : 
Tinker et al.ll2008l : iMore et al.il201ll '). 



In our analysis we apply the SO method, with masses 
measured at four different overdensities corresponding to 
Ac = 200, 500 and 1500, thus ranging from the overden- 
sity which characterize the whole virialized region of halos 
up to the typical overdensity which is accessible by Chandra 
and XMM-Newton X-ray observations. Our halo identifica- 
tion proceeds in two steps. In the first step we run a FoF 
algorithm with linking length 6 = 0.16 over the distribu- 
tion of DM particles (in the DM-only simulation the FoF is 
run over the distribution of more massive particles). Then, 
we identify in each FoF group the DM particle which cor- 
responds to the minimum of the potential. The position of 
this particle is taken to be the center of the cluster from 
where to grow spheres whose radius is increased until the 
mean density within it reaches the required overdensity Ac. 
The mass Ma^ within this spherical region of radius i?Ac is 



4 3 



(1) 



Since each halo is firstly identified starting from a FoF 
algorithm, it inherits some FoF disadvantages. A well known 
potential problem with FoF is that there are situations in 
which two halos are connected through a bridge of particles. 
Since this halo is counted only once, this could affect the 
number of SO halo number and the resulting mass function. 



As discussed by jRee d et al.ll2007l ) , this effect becomes more 
important at high redshift and for poorly resolved low-mass 
halos. Since we restrict our analysis to halos having a min- 
imum mass of 10^"^ Mq, thus being resolved by at least 
10^ particles) and redshift 2; ^ 1, and we use a linking length 
smaller than the usually adopted value b — 0.2, we expect 
the bias induced by using FoF parent groups should be miti- 
gated. Since the FoF grouping is carried out using DM parti- 
cles as primary particles, we expect halo bridging to affect in 
the same way the N-body and the different hydrodynamical 
simulations. Therefore, our main conclusions on the rela- 
tive effect of baryons on the mass function should be left 
unchanged by the effect of using FoF groups as the starting 
point of the SO identification. Finally, since the groups iden- 
tified by FoF algorithm have by definition no overlapping, 
we do not include in our identification o f SO halos any re- 
striction to prevent such overlapping (see [Tinker et al. H2OO8I 
for a discussion on halo overlapping). 

3.2 Effect on halo mass and density profile 

We first focus on the impact that baryons have on the mass 
of individual halos. To this purpose, we show in Figure [2] 
the distribution of the differences between halos identified 
in the two hydrodynamical simulations and in the DM sim- 
ulation, at different overdensities. Results in this figure are 
shown for all the halos that in the DM simulation have 
Ma ^ 10^'^ Mq. To compare halos masses, one has to 
identify a halo selected in the DM simulation with its coun- 
terpart in each one of the hydrodynamical simulations. The 
easiest way to perform this identification is to look for the 
halos having the closest coordinates. While this procedure 
provides a reliable identification of corresponding halos in 
two different simulations for the most massive systems, it 
turns out not to be accurate for poorer systems. In fact, 
besides affecting the mass of halos, the presence of baryons 
also slightly alter the overall dynamics and, therefore, the 
exact halo positions. In order to overcome this difficulty we 
decided to follow a different procedure to find in each of 
the hydrodynamical simulations the halos corresponding to 
those identified in the DM run. For each halo in the DM 
simulation we identify the Lagrangian region from where 
particles following within its virial radius by 2 = come 
from. We then look in each of the GH and CSF simulations 
for a halo that contains at least 60 per cent of the particles 
coming from the same Lagrangian region. We verified that 
the final results do not change significantly if we use instead 
a more restrictive requirement to find instead 80 per cent of 
the particles from the same Lagrangian region. 

From Fig. [2l we see that significant mass differences, of 
up to 20 per cent, are found for Ac — 1500, with the distri- 
bution of such differences becoming narrower at Ac — 200. 
Correspondingly, the mean value of the halo mass increase 
induced by the presence of baryons decreases from ~ 6-7 
per cent at Ac = 1500 to ~ 3-4 per cent at Ac — 500, while 
being ^ 1 per cent at Ac — 200. Furthermore, any differ- 
ences between the two hydrodynamic runs is much smaller 
than the difference that each of them has with respect to 
the DM run. This result is in line with the weak sensitivity 
of the mass function on the details of the baryon physics, as 
shown in Fig. [5] Even at the highest considered overdensity, 
Ac — 1500, there is a small number of halos whose mass in 
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Figure 2. The distribution of the mass ratios between the hales identified in each of the two hydrodynamical simulations to the 
corresponding halos in the DM simulation. The three panels from left to right show the results at Ac = 1500, 500, 200. Red and black 
continuous histograms show the results for Mqu /Mqm and Mcgp /Mj^m y respectively. In each panel, results are shown only for halos 
with ^ 10^*/i~^ M0. The solid vertical lines correspond to no mass variation, while the dashed and dotted vertical lines show the 

mean values (which is shown on the right-top of each panel) of Mqu /M]ji^j and McsF / ^DM ) respectively. 



the hydrodynamical runs is smaller than in the DM run. The 
reason for this is the different timing of merging of substruc- 
tures in the different simulations. This occasionally causes 
some of these substructures to be found outside -Risoo while 
located within this radius in the DM simulation. 

In order to quantify a possible mass dependence of this 
halo mass difference, we show in Figure [3] the mean value 
of such a difference for each mass bin where the mass func- 
tion is computed. Throughout our analysis, we use a fixed 
mass bin with width AlogAf = 0.2. With such a narrow 
bin, the mass function suffers for large sampling effect in 
the high-mass end, due to the exponential dearth of the 
massive halo population. To overcome such sampling effect, 
we merge mass bins containing less than 10 objects into the 
adjacent lower mass bin. Each mass bin is then weighted 
proportionally to the number of clusters it contains. 

The increase of halo masses in both the GH and CSF 
simulations is to good approximation independent of halo 
mass, at least for log(M//i~^ M©)^ 13.5, at overdensities 
Ac = 200 and 500. Again, this shift in mass turns out to be 
similar in the two hydrodynamical runs. It amounts to about 
1-2 per cent at Ac = 200 and ~ 4 per cent at Ac = 500, 
in line with the results shown in Fig. (2] The increasing star 
formation efficiency in lower mass halos makes the mass in- 
crease in the CSF simulation to be larger than for the GH 
case. This difference between GH and CSF halo masses fur- 
ther increases for Ac = 1500. At this overdensity we can not 
define a mass range over which the increase of halo masses 
due to baryons is nearly constant and independent of gas 
physics. 

To better understand the origin of the mass difference 
between halos identified in different simulations, we further 
show in Figure |4] the radial profile of the mean total den- 
sity for halos identified in the three simulations. The four 
panels correspond to different mass ranges. Since density 
is normalized to p200, i-e. the mean density within -R200, 
the profiles reach the unity value for R/R200 = 1. As for 
the halos identified in the GH simulation (red dot-dashed 
curves), their profiles have small but sizable differences with 



respect to the DM case (solid black curves). At intermedi- 
ate radii, 0.1^ -R/ii2oo~ 1 the GH profiles lie above those of 
the DM simulation. This result, which holds in dependent of 
the ha lo mass, is consistent with that found bv lRasia et al.l 
(2004|) in their comparison of halo profiles from DM-only 
and non-radiative hydrodynamical simulations. These au- 
thors argued that the more concentrated density profiles in 
non-radiative simulations, with respect to DM-only simu- 
lations, is the result of energy redistribution between the 
DM and the baryo nic component during halo collapse (see 
also iLin et akllifTO ) . We postpone to a forthcoming paper a 
detailed comparison between concentrations for halos iden- 
tified in DM and hydrodynamic simulations (Rasia et al. in 
preparation). It is only at small radii, R ^ O.O8-R200, that 
gas pressure support makes the total density profiles in the 
GH simulation slightly flattening with respect to the DM 
simulation. 

As for the radiative CSF simulation (blue dashed 
curves), the sinking of cooled baryons, converted into stars, 
in the central halo regions causes the already known effect of 
adiabatic contraction, with the resulting steepening of the 
density profiles in these regions. A comparison of the re- 
sulting profiles for the different mass ranges indicates that 
this effect is more pronounced for halos of smaller mass, 
consistent with the expectation that cooling is in fact more 
efficient in lower mass halos, due to their higher concen- 
tration. The effect of gas cooling is rather pronounced for 
R^ 0.2_R2oo- We note that the vertical purple line in Fig- 
ure |4] mark the mean value of -R1500 for halos of the DM 
simulation. The corresponding value of -Risoo for the CSF 
simulation is in fact slight larger than in the DM case. This 
difference explains why halo masses in the CSF simulation 
are only slightly larger than in the GH simulation already 
at -R1500. 



3.3 Effect on the Halo Mass Function 

In order to compute the mass function, we group SO halos 
within mass bins having fixed width A log M = 0.2. Then, 



© 0000 RAS, MNRAS 000, 000-000 



6 Weiguang Cui, et al. 



1 .20 



1.15 



1.10 



1 .05 



.00 



A =500 



A^ = 200 




fc-&-^fi~"n — i — ^ — ^ H — ^ 



12.5 13.0 13.5 14.0 14.512.5 13.0 13.5 14.0 14.5 

log(M,, /h-^Me) log(M,, /K^U^) 



Figure 3. The ratio between masses of halos identified in tfie hydrodynamical simulations and tfic corresponding halos from the DM 
simulation, as a function of halo mass in the DM run, Mj^m- Left and right panels show the results for the GH and CSF run, respectively. 
In each panel different line styles, associated with different colors, corresponds to different redshifts. Different symbols indicate instead 
different Ac. For reference, the horizontal light long-dashed line correspond to no mass variation. 



the mass assigned to each bin is computed as the mean over 
all the halos belonging to that mass bin. Whatever proce- 
dure one adopts to choose the mass to be assigned to a 
given bin, it is clear that the binning procedure introduce 
an uncertainty i n the resulting mass function. As discussed 
bv lLukic et all ^Mh . this uncertainty is negligible as log 
as the bin width does not exceed AlogM = 0.5. 

We show in Figure [5] the HMF for our three simula- 
tions, computed for Ac = 200, 500, 1500 (from upper to 
lower groups of curves). To better emphasize the mass vari- 
ation induced by the presence of baryons, we show in the 
lower panel the difference in the number of clusters within 
each mass bin, between each of the two hydrodynamical sim- 
ulations and DM simulation. Due to the specific treatment 
of the last bin, its width can be different for different simu- 
lations. Therefore, when we compare the number of objects 
in such last bins, we rescale the cluster counts within each of 
them by scaling it to the bin width in units of A log M — 0.2. 

In general, we find that the presence of baryons leads to 
an increase of the HMF by an amount increasing as we move 
to more internal regions at higher Ac. In general, this vari- 
ation is nearly independent of mass, except possibly in the 
high mass end, beyond log(M/^~^ M©) ~ 14.5. This is the 
regime where exponential tail takes place. Given the limited 
box size, the resulting limited statistics of massive halos does 
not allow us to draw robust conclusions for such high masses, 
especially when considering Ac — 500 and above. For the 
GH non-radiative simulation the HMF increase is negligible 
at Ac = 200, and amounts to ^ 3 per cent at the largest 
sampled masses. This difference increases to ^ 8 per cent as 
we move to Ac = 500, at least up to log{M/h~^ Mq) ~ 14.5. 
We note that at such overdensities the effect of introducing 
baryons produce a variation of the HMF with respect to the 



DM simulation which is larger than the difference between 
the GH and the CSF run. This indicates that, while it is im- 
portant to account for the presence of baryons in the HMF 
calibration for Aci$ 500, the details of the physical processes 
regulating their evolution has a minor impact. At a higher 
overdensity Ac = 1500, the effect of radiative physics is of 
increasing the HMF by about 20 per cent for CSF run and 
around 10 per cent for GH run. This result is in line with the 
expectation that a more concentrated density profile in the 
presence of gas cooling jCnedin et al. l2004l :[ Pedrosa et al.l 
l2009l : iTissera et ai]|2010l : iDuffv et al-lbOld ) 

In general, our results for an increase of the HMF for the 
GH simulation is in line with previous findings, also based 
on non-radiative simulations, for an inc rease of halo con- 
centration induced by the presence of gas ([Rasia et al.ll2004l : 
iJing et all [200^ : iLin et al.|[200^ : iRudd et al.ll2008l ). The ef- 
fect of the physics of baryons becomes more important at 
Ac = 1500. I n the ir analysis of the cumulative mass function 
iRudd et al.l (|2008l ') found that the presence of non-radiative 
gas induces a negligible HMF variation for masses estimated 
at the virial radius, corresponding to Ac — 100 for their sim- 
ulated cosm o logy. While this result is in agreement with our, 
iRudd et all (|2008l 'l find that the HMF increases by about 
10 per cent when radiative cooling and star formation are 
included. One possible reason for the different effect of ra- 
diativc physics in our analysis and in that of iRudd et al.l 
(2008) could lie in the different efficiency of the feedback 
included in the simulations. In our case we include a rather 
efficient SN feedb ack, that could mitig ate the effects of adia- 
batic contraction. Istanek et al.l (|2009l ) compared results for 
a non-radiative simulation and for a pre-heated radiative 
simulations. They found that at A = 500 the latter predicts 
a HMF which is lower than the former. The reason for this 
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Figure 4. Mean radial profiles of the density within a given radius. Density is expressed in units of P20O1 which is the mean density within 
R200- In each panel, solid (black), dot-dashed (red) and dashed (blue) curves correspond to the DM, GH and CSF runs, respectively. 
The two vertical lines mark the mean value of -R1500 and ^500 for the DM simulation (purple triple-dot-dashed and green long-dashed 
lines, respectively). The four panels correspond to four different mass ranges over which the mean profiles are computed, as indicated in 
the labels. 



is that the fairly strong pre-heating introduced in their sim- 
ulation at z = 4 devoid halos by a substantial amount of 
gas, which is later prevented to re-accelerate in the forming 
larger halos. 

iTinker" et al.l (|2008l ) investigated the redshift evolution 
of the mass function computed at different values of Ac 
based on simulations including only dark matter. They 
found that the abundance of halos at a given logcr"^ mono- 
tonically decreases with increasing z in the interval [0,2.5], 
where a is r.m.s. variance of the linear density field smoothed 
on a given mass scale. In the following we will discuss the 
impact that baryons in the GH and CSF simulations have 
on the evolution of the HMF. 

We show in Figure[6]the evolution of the HMF at three 
different redshifts, z = 0, 0.6 and 1, for Ac = 500. Similarly 
to Fig. [S] we also show in the bottom panel the ratio between 
the HMF from each of the two hydrodynamical simulations 
and that of the DM simulation. In general, we find that 
the effect of baryons on the MF is to good approximation 
independent of mass at all redshifts for Ac — 500. As for 
the non-radiative GH simulation, the increase of the HMF 
with respect to the DM case is always very small out to 
z — 1, and ;$ 8 per cent. A more significant effect is instead 
found for the radiative CSF simulation. In this case the HMF 



Redshift 


GH 




CSF 






Ac = 500 


Ac = 200 


Ac = 500 


Ac = 200 


z = 0.0 


1.037 


1.012 


1.049 


1.016 


z = 0.6 


1.028 


1.000 


1.044 


1.005 


2 = 1.0 


1.028 


0.997 


1.040 


0.999 



Table 1. Mean values of the ratio between halo masses in the 
hydrodynamical and in DM simulations, at different redshifts and 
Ac values, for both the non-radiative (GH) and radiative (CSF) 
simulations. Such values have been computed by including only 
halos withlog(MA„//i~^ Mq) ^ 13.5. 



increases by a larger amount at progressively higher redshift, 
reaching the ~ 10 per cent level at 2; = 1. This result agrees 
with the expectation that a more efficient cooling takes place 
at higher redshift, which induces a stronger effect on halo 
masses. 

Due to the complexity and partial knowledge of the 
baryonic process taking place in galaxy formation, accu- 
rately calibrating the HMF from hydrodynamical simula- 
tion may appear a challenging task. However, owing to the 
results show in Fig. [S] it turns out that variation of halo 
masses in both GH and CSF simulations, with respect to 
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Figure 7. The number difference after iialo mass calibration. Different colorful lines show the redshift, while the two line styles, solid 
and dotted represent Ac = 500, 200, respectively. The left panel shows the corrected results for GH run, and the right is for CSF run. 



the pure DM case, is to good approximation independent of 
halo mass, at least for log(M//i"^ M©) ^ 13.5 and Ac = 200 
and 500. Therefore, we can in principle attempt to account 
for the effect of baryons in the HMF by a simple shift of 
halo masses, at least within the mass range and for the Ac 
values for which is approximation is expected to hold. In 
Table [T] we report the average values of the mass ratio for 
halos identified in the hydrodynamical and in the DM sim- 
ulations, for different Ac and redshift values. This values, 
which are computed as an average over all halos having mass 
log(MA,//i-'M0) ^ 13.5, are then used to rescale the HMF 
from the hydrodynamical simulation. 

In Figure [T] we show the ratio between the number of 
halos of different mass found in the hydrodynamical and 
in the DM simulations, after applying the mass shifts re- 
ported in Table [1] A larger mass binning is used here, 
AlogM = 0.3, to reduce fluctuations associate to sampling 
noise. From the left panel of Figure [71 the difference in the 
halo number is now consistent with zero, with fluctuations 
around this value of ^ 3 percent for Ac — 500. This result 
holds independently of mass and redshift, at least for halos 
with log(A/A^//i~^ Mq) 13.5. As expected, the correc- 
tion is less effective for smaller masses, owing to the larger 
mass difference induced by baryonic effects in smaller halos. 
Clearly, the correction is less pronounced at Ac = 200, owing 
to the smaller impact of baryons at this overdensity. How- 
ever, also in this case, correcting the HMF according to a 
unique mass shift further reduces the difference between hy- 
drodynamical and DM simulations at the 1-2 per cent level. 
The number difference for CSF run at Ac = 500 is also sup- 
pressed to unity for all halos with \og{M /h~^ M©) ^ 13.5. 

In conclusion, the results obtained from our analysis 
indicates that the relative variation of halo masses due to 
baryon effects are always within 5 per cent, for both non- 
radiative and radiative simulations, also almost independent 
of redshift. This result holds for masses computed at over- 
density Ac = 200 and 500, and for halos having mass at 
least comparable to that of a galaxy group. Correcting the 



mass function with a constant mass shift in this mass range 
largely accounts for the differences between hydrodynamical 
and DM simulations. 



4 DISCUSSION AND CONCLUSION 

In this paper we presented an analysis of the effect of baryons 
on the calibration of the halo mass function (HMF). To this 
purpose, we carried out one DM-only simulation (DM) and 
two hydrodynamical simulations, a non-radiative one in- 
cluding only the effect of gravitational gas heating (GH) and 
a radiative one including also the effect of star formation and 
SN feedback in the form of galactic ejecta. The three simula- 
tions, which_axejll_based on the Tree-PM/SPH GADGET-3 
code (|Springej|2005al '). started from exactly the same initial 
conditions and followed the evolution of 2 x 1024'' particles 
within a box having a comoving size of 410 h-^ Mpc. Ha- 
los have been identified using a spherical overdensity (SO) 
algorithm, and results have been presented at three red- 
shifts, 2 = 0, 0.6 and 1. Halo masses have been computed 
at different overdensities (with respect to the critical one), 
Ac = 200, 500 and 1500. The main results of our analysis 
can be summarized as follows. 

1. The fractional difference between halo masses in the 
hydrodynamical and in the DM simulations is found to 
be almost constant, at least for halos more massive than 
log(MAc//i~^ M0) ^ 13.5. In this range, the mass increase 
in the hydrodynamical simulations is of about 4-5 per cent 
at Ac = 500 and 1-2 per cent at Ac = 200. Quite interest- 
ingly, these differences are nearly the same for the GH and 
the CSF simulations (see Fig|3]and Table [T|. Such relative 
mass variations can not be considered any more as constant 
at higher overdensity, Ac = 1500, and smaller masses. In 
these cases, mass difference markedly increases for smaller 
halos in the CSF simulation, while it decreases in the non- 
radiative GH simulation. 

2. These variations of halo masses induce corresponding 
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Figure 5. The halo mass function for our simulations, with 
masses computed at different overdensities Ac. Results for the 
DM simulations are shown with dotted curves, while results for 
the GH and CSF hydrodynamical simulations are shown with 
the dashed and dot— dashed curves, respectively. Upper to lower 
curves correspond to result for Ac = 200, 500 and 1500 (green, 
red and black curves, respectively). The lower panel shows the 
ratio between the number of halos found in each mass bin for 
each of the two hydrodynamical simulation and the DM simu- 
lation. We apply a linear interpolation of the mass functions to 
compute the difference in the halo number exactly for the same 
mass values. 



variations of the HMF (see Fig. [5}. At « = 0, the HMFs 
for GH and CSF simulations are close to the DM one, with 
differences of ^ 3 per cent at Ac — 200, in line with the 
small correction in halo masses. Such a difference increases 
to ~ 7 per cent at Ac = 500 and reaches ^ 10-20 per 
cent at Ac = 1500. At the latter overdensity, the increase 
in the HMF for the CSF run is larger by about a factor 
2 with respect to the GH run. This result in line with the 
expectation that baryonic processes have a stronger impact 
in the central halo regions. At higher redshift, differences 
with respect to the DM HMF tend to increase, especially 
for Ac = 1500 (see [S| and for the CSF case. Again, this 
result agrees with the increase of cooling efficiency within 
halos at higher redshift. 

3. Based on the above results, we showed that assuming a 
constant mass variation to rescale the HMF from the hydro- 
dynamic simulations reduces the difference with respect to 
the DM case. We apply a uniform mass shift, calibrated for 
halo masses log(M fe"^ Mq) ^ 13.5 for Ac = 200 and 500. 
We verified that the difference between hydrodynamical and 
DM HMFs becomes negligible, with fluctuations around null 
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Figure 6. The redshift evolution of the halo mass function for 
our simulations, with masses computed overdensity Ac = 500. 
Results for the DM simulations are shown with dotted curves, 
while results for the GH and CSF hydrodynamical simulations 
are shown with the dashed and dot-dashed curves, respectively. 
Upper to lower curves correspond to result for z = 0.00, 0.58 and 
1.00 (black, red and green curves, respectively). The lower panel 
shows the ratio between the number of halos found in each mass 
bin for each of the two hydrodynamical simulation and the DM 
simulation. 



of ^ 3 per cent at Ac = 500. Even though mass variations 
are smaller at Ac = 200, we still find that a uniform mass 
rescaling gives a small but sizable reduction of the HMF 
difference also at this overdensity. 

The future generation of large surveys of galaxy clus- 
ters, from X-ray, optical and Sunyaev-Zeldovich observa- 
tions, could provide stringent constraints of cosmological 
parameters through the study of the evolution of the mass 
function. However, a necessary condition to fully exploit the 
cosmological information content of such surveys is that the 
theoretical mass function needs to be calibrated to a preci- 
sion better than 10 per cent (e.g. IWu et al]|2010l ). In this 
respect, the results of our analysis have interesting implica- 
tions to gauge the uncertainty in the mass function calibra- 
tion associated to the uncertain baryon physics. 

First of all, the HMF turns out to be less prone to 
such effects if computed at Ac = 500, while they become 
more important and likely difficult to model in detail at 
higher overdensities. Furthermore, adopting a constant mass 
shift provides a rather accurate correction to the HMF cali- 
brated from DM simulations, at least for halos having size or 
galaxy groups or larger. This result holds for both the non- 
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radiative (GH) and the radiative (CSF) simulations, which 
have rather similar mass corrections at Ac = 500. Since 
the CSF run only include SN feedback, but no AGN feed- 
back, clusters in this simulations still suffers for overcooling. 
Therefore, the GH and CSF simulations should in princi- 
ple bracket the case in which the correct amou nt of baryons 
cools within DM halos. However, we note that lStanek et al] 
found a slight decrease, rather than an increase of 
the HMF in simulations including an impulsive pre-heating. 
Since a phenomenological pre-heating only provides an ap- 
proximate description of the astrophysical mechanisms reg- 
ulating star formation, it would be interesting to repeat our 
analysis also in the presence of a mechanism for AGN feed- 
back that regulates cooling in groups and clusters to the ob- 
served level (e.g ., jPuchwcin ct al. 2008; Fabian et al. 20101: 
iMcCarthv et al.i[20iq ). 

Another direction in which our analysis should be im- 
proved concerns the size of simulation box, so as to better 
sample the population of massive halos. Although our re- 
sults indicate that a mass-independent mass shift should 
be applied to account for baryonic effects, one may won- 
der whether this prescription can be extrapolated to the 
most massive halos, whose population is mostly sensitive 
to choice of the cosmological model. Future development in 
supercomputing capabilities will soon open the possibility 
to carry out hydrodynamical simulations which will cover 
dynamic ranges comparable to those accessible by the N- 
body simulations currently used to calibrate the halo mass 
function. 
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